#include "AnaDef.h"
#include "XGLUtils.h"

const Double_t EPSILON=1e-12;

void renormalizeWithoutE(const TString  fin)
{
  TFile *ff = new TFile(fin);

  TGraphAsymmErrors * gstat[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * gsys[]={0x0,0x0,0x0,0x0};
  const TGraphAsymmErrors * gdumm[]={0x0,0x0,0x0,0x0};
  for(Int_t itype=0; itype<4; itype++){
    gstat[itype]=(TGraphAsymmErrors *)ff->Get(Form("grstat%d",itype));
    gsys[itype]=(TGraphAsymmErrors *)ff->Get(Form("grsys%d",itype));
    if(!gstat[itype] || !gsys[itype]){
      printf("no stat or sys %d\n", itype);exit(1);

    }

    gstat[itype]->SetTitle("Normalized to sum (pi, k, p, e)");
    gsys[itype]->SetTitle("Normalized to sum (pi, k, p, e)");

    gdumm[itype]=(TGraphAsymmErrors *)gstat[itype]->Clone(Form("dum%d",itype));
  }

  const Int_t ngr=4;
  Int_t id0s[ngr];
  Int_t id1s[ngr];
  const Int_t nptmp = XGLUtils::GetCommonLimits(ngr, gdumm, id0s, id1s);

  /*
  const Int_t np=gstat[0]->GetN();
  if(np!=gstat[1]->GetN() || np!=gstat[2]->GetN() || np!=gstat[3]->GetN()){
    printf("bad np!! %d %d %d %d\n", np, gstat[1]->GetN(), gstat[2]->GetN(), gstat[3]->GetN()); exit(1);
  }
  */

  const Int_t np = nptmp;
  for(Int_t rawip=0; rawip<np; rawip++){
    Int_t ips[4];
    for(Int_t it=0; it<4; it++){
      ips[it] = id0s[it] + rawip;
    }

    const Double_t xx = gstat[0]->GetX()[ips[0]];
    if(fabs(xx-gstat[1]->GetX()[ips[1]])>EPSILON){
      printf("bad xx %d %e %e\n", rawip, xx, gstat[1]->GetX()[ips[1]]); exit(1);
    }

    //0, 1, 2 only proton, pion and kaon
    Double_t sum3 = 0;
    for(Int_t it=0; it<3; it++){
      sum3 += gstat[it]->GetY()[ips[it]];
    }

    //0, 1, 2, 3 including electron
    const Double_t sum4 = sum3 + gstat[3]->GetY()[ips[3]];

    const Double_t corr = sum4/(sum3+EPSILON);

    printf("ipoint %d sum3 %e sum4 %e corr %e before %e %e %e %e\n", rawip, sum3, sum4, corr, gstat[0]->GetY()[ips[0]], gstat[1]->GetY()[ips[1]], gstat[2]->GetY()[ips[2]], gstat[3]->GetY()[ips[3]]);

    //only for itype 0, 1, 2
    for(Int_t itype=0; itype<=2; itype++){
      gstat[itype]->GetY()[ips[itype]] *= corr;
      gstat[itype]->GetEYlow()[ips[itype]] *= corr;
      gstat[itype]->GetEYhigh()[ips[itype]] *= corr;

      gsys[itype]->GetY()[ips[itype]] *= corr;
      gsys[itype]->GetEYlow()[ips[itype]] *= corr;
      gsys[itype]->GetEYhigh()[ips[itype]] *= corr;

      gstat[itype]->SetTitle("Renormalized to sum (pi, k, p)");
      gsys[itype]->SetTitle("Renormalized to sum (pi, k, p)");
    }

    printf("ipoint %d sum3 %e sum4 %e corr %e after  %e %e %e %e\n", rawip, sum3, sum4, corr, gstat[0]->GetY()[ips[0]], gstat[1]->GetY()[ips[1]], gstat[2]->GetY()[ips[2]], gstat[3]->GetY()[ips[3]]);
  }

  TString fname(fin);
  while(fname.Contains("/")){
    fname=fname(fname.First('/')+1,fname.Length());
  }

  fname.Prepend("Renormalized_");

  TFile *fout = new TFile(fname,"recreate");
  for(Int_t itype=0; itype<4; itype++){
    gstat[itype]->Write();
    gsys[itype]->Write();
  }
  fout->Save();
  fout->Close();
  delete fout;
}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=2!!");return 1;
  }

  renormalizeWithoutE(argv[1]);
  return 0;
}
